import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as snsNa początek ustawimy globalnie opcję set_output ze
Scikit-learn'a (dokumentacja).
Dzięki temu transformowane dane będą zawsze na wyjściu typu
pd.DataFrame, zamiast np.ndarray, i zachowają
odpowiednie nazwy cech.
import sklearn
sklearn.set_config(transform_output="pandas")Wykorzystamy zbiór danych Ames housing, w którym zadaniem jest przewidywanie wartości domu na podstawie cech budynku, działki, lokalizacji itp. Jest to więc przewidywanie wartości ciągłej, czyli regresja. Obejmuje zmienne do usunięcia (np. ID transakcji), numeryczne (floaty i inty), kategoryczne nieuporządkowane (categorical nominal) oraz kategoryczne uporządkowane (categorical ordinal), więc będzie wymagał wstępnego przetworzenia tak jak większość prawdziwych danych w ML.
Inne znane, ale gorsze jakościowo zbiory tego typu to na przykład:
Autor zbioru to Dean De Cock, a zbiór został opisany oryginalnie w tym artykule.
Szczegółowe opisy zmiennych znajdują się w pliku
ames_description.txt.
Formatem pliku jest Apache Parquet, bardzo wydajny format binarny. Ma on szereg zalet w ML i inżynierii danych porównaniu do plików CSV:
Dla zainteresowanych szczegółami: link 1, link 2, link 3. Ten format ma szerokie wsparcie, w tym w Pandasie. W przypadku tego zbioru, CSV zajmuje ~1.2 MB, a Parquet ~0.2 MB.
df = pd.read_parquet("ames_data.parquet")
# remove dots from names to match data_description.txt
df.columns = [col.replace(".", "") for col in df.columns]
df.head()| Order | PID | MSSubClass | MSZoning | LotFrontage | LotArea | Street | Alley | LotShape | LandContour | ... | PoolArea | PoolQC | Fence | MiscFeature | MiscVal | MoSold | YrSold | SaleType | SaleCondition | SalePrice | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | 526301100 | 20 | RL | 141.0 | 31770 | Pave | None | IR1 | Lvl | ... | 0 | None | None | None | 0 | 5 | 2010 | WD | Normal | 215000 |
| 1 | 2 | 526350040 | 20 | RH | 80.0 | 11622 | Pave | None | Reg | Lvl | ... | 0 | None | MnPrv | None | 0 | 6 | 2010 | WD | Normal | 105000 |
| 2 | 3 | 526351010 | 20 | RL | 81.0 | 14267 | Pave | None | IR1 | Lvl | ... | 0 | None | None | Gar2 | 12500 | 6 | 2010 | WD | Normal | 172000 |
| 3 | 4 | 526353030 | 20 | RL | 93.0 | 11160 | Pave | None | Reg | Lvl | ... | 0 | None | None | None | 0 | 4 | 2010 | WD | Normal | 244000 |
| 4 | 5 | 527105010 | 60 | RL | 74.0 | 13830 | Pave | None | IR1 | Lvl | ... | 0 | None | MnPrv | None | 0 | 3 | 2010 | WD | Normal | 189900 |
5 rows × 82 columns
df.shape(2930, 82)
Wstępne czyszczenie danych (data cleaning) obejmuje:
Order oraz PIDGrnHill oraz
Landmrk, które obejmują w sumie zaledwie 3 domyTo drugie jest motywowane wykresem przedstawionym poniżej. Zostało uznane za błąd już przez samego autora.
df = df.drop(["Order", "PID"], axis="columns")
df = df.loc[~df["Neighborhood"].isin(["GrnHill", "Landmrk"]), :]
plt.scatter(df["GrLivArea"], df["SalePrice"])
plt.title("House area vs price")
plt.xlabel("GrLivArea")
plt.ylabel("SalePrice")
plt.show()
df = df.loc[df["GrLivArea"] <= 4000, :]
plt.scatter(df["GrLivArea"], df["SalePrice"])
plt.title("House area vs price, outliers removed")
plt.xlabel("GrLivArea")
plt.ylabel("SalePrice")
plt.show()

Zawsze warto też przyjrzeć się rozkładowi zmiennej docelowej, żeby poznać jego typ i skalę. Jak widać poniżej, rozkład jest dość skośny, co ma sens - mało jest bardzo drogich domów. Można to zatem uznać za problem regresji niezbalansowanej (imbalanced regression), w którym rozkład zmiennej docelowej nie jest normalny.
Załóżmy, że typowego klienta nie interesują bardzo drogie domy. Cały ciężki ogon rozkładu można zatem uznać zasadniczo za outliery - im większa wartość, tym mniej interesuje przeciętną osobę. Warto też pamiętać, że realny "przeciętny" dochód to mediana, a nie średnia zarobków, co odpowiednio przekłada się na zainteresowanie domami.
W tym przypadku warto przetransformować tę zmienną logarytmem bliżej rozkładu normalnego z dwóch powodów.
Po pierwsze, stabilność numeryczna - wartości są bardzo duże, rzędu setek tysięcy, i przy obliczaniu odchylenia standardowego czy wariancji to byłby problem. Logarytm to przykład transformacji stabilizującej wariancję (variance-stabilizing transform).
Ponadto regresja liniowa (pierwszy model liniowy, który będziemy rozważać) czyni szereg założeń w modelowaniu, co wynika z jej podstaw statystycznych. Są to między innymi:
W naszym przypadku oznacza to, że podobnie często nie doceniamy albo przeceniamy wartość domu, oraz dla tanich i drogich domów mamy podobny błąd. Przy rozkładzie zmiennej zależnej bardzo różnym od normalnego te założenia są typowo łamane, co skutkuje gorszą jakością predykcji modelu. Dlatego transformacji warto też dokonać dla lepszej jakości modelu.
Skąd wiemy, że logarytm? Bo rozkłady o dodaniej skośności (positively skewed) przesuwa się "w prawo" właśnie logarytmem, a w drugą stronę eksponentą. Dla zainteresowanych: link 1, link 2. Są też algorytmy, które potrafią same estymować potrzebną transformację, np. Box-Cox transform albo Yeo-Johnson transform. Podobne transformacje można też oczywiście wykonywać nie tylko dla zmiennej zależnej, ale też (a nawet jest to częstsze) dla cech.
Zadanie 1 (1 punkt)
"SalePrice".
Pamiętaj o podaniu tytułu i dodaniu legendy.np.log1p (zamień istniejące
wartości).Może ci się przydać metoda .plot.hist() dla
pd.Series. Zwraca ona normalny wykres Matplotliba, do
którego można potem dodawać kolejne elementy.
# Plot before taking logarithm
plt.figure(figsize=(10, 6))
df["SalePrice"].plot.hist(edgecolor='black', bins=30, alpha=0.7, label="Sale Price")
plt.title("Sale Price Histogram")
plt.xlabel("Sale Price")
plt.ylabel("Number of holders")
plt.legend()
mean_price = np.mean(df["SalePrice"])
median_price = np.median(df["SalePrice"])
print(f"Średnia: {int(mean_price)}")
print(f"Mediana: {int(median_price)}")
plt.axvline(mean_price, color='red', linestyle="--")
plt.axvline(median_price, color='blue', linestyle="--")
plt.text(mean_price, 0 , "mean", color="red", rotation=60, fontsize = 12)
plt.text(median_price, 100 , "median", color="blue", rotation=60, fontsize = 12)
# Plot after taking logarithm
df["SalePrice"] = np.log1p(df["SalePrice"])
plt.figure(figsize=(10, 6))
df["SalePrice"].plot.hist(edgecolor='black', bins=30, alpha=0.7, label="Sale Price")
plt.title("Sale Price Histogram after taking logarithm")
plt.xlabel("Sale Price")
plt.ylabel("Number of holders")
plt.legend()
mean_price_log = np.mean(df["SalePrice"])
median_price_log = np.median(df["SalePrice"])
print(f"Średnia (po zlogarytmowaniu): {int(mean_price_log)}")
print(f"Mediana (po zlogarytmowaniu): {int(median_price_log)}")
plt.axvline(mean_price_log, color='red', linestyle='--')
plt.axvline(median_price_log, color='blue', linestyle='--')
plt.text(mean_price_log, 0 , "mean", color="red", rotation=60, fontsize = 12)
plt.text(median_price_log, 100 , "median", color="blue", rotation=60, fontsize = 12)
plt.show()Średnia: 180358
Mediana: 160000
Średnia (po zlogarytmowaniu): 12
Mediana (po zlogarytmowaniu): 11


W zbiorze znajdują się liczne wartości brakujące, które zostały uzupełnione na podstawie tego notebooka na Kaggle.
def replace_na(df: pd.DataFrame, col: str, value) -> None:
df.loc[:, col] = df.loc[:, col].fillna(value)
# Alley : data description says NA means "no alley access"
replace_na(df, "Alley", value="None")
# BedroomAbvGr : NA most likely means 0
replace_na(df, "BedroomAbvGr", value=0)
# BsmtQual etc : data description says NA for basement features is "no basement"
replace_na(df, "BsmtQual", value="No")
replace_na(df, "BsmtCond", value="No")
replace_na(df, "BsmtExposure", value="No")
replace_na(df, "BsmtFinType1", value="No")
replace_na(df, "BsmtFinType2", value="No")
replace_na(df, "BsmtFullBath", value=0)
replace_na(df, "BsmtHalfBath", value=0)
replace_na(df, "BsmtUnfSF", value=0)
# Condition : NA most likely means Normal
replace_na(df, "Condition1", value="Norm")
replace_na(df, "Condition2", value="Norm")
# External stuff : NA most likely means average
replace_na(df, "ExterCond", value="TA")
replace_na(df, "ExterQual", value="TA")
# Fence : data description says NA means "no fence"
replace_na(df, "Fence", value="No")
# Functional : data description says NA means typical
replace_na(df, "Functional", value="Typ")
# GarageType etc : data description says NA for garage features is "no garage"
replace_na(df, "GarageType", value="No")
replace_na(df, "GarageFinish", value="No")
replace_na(df, "GarageQual", value="No")
replace_na(df, "GarageCond", value="No")
replace_na(df, "GarageArea", value=0)
replace_na(df, "GarageCars", value=0)
# HalfBath : NA most likely means no half baths above grade
replace_na(df, "HalfBath", value=0)
# HeatingQC : NA most likely means typical
replace_na(df, "HeatingQC", value="Ta")
# KitchenAbvGr : NA most likely means 0
replace_na(df, "KitchenAbvGr", value=0)
# KitchenQual : NA most likely means typical
replace_na(df, "KitchenQual", value="TA")
# LotFrontage : NA most likely means no lot frontage
replace_na(df, "LotFrontage", value=0)
# LotShape : NA most likely means regular
replace_na(df, "LotShape", value="Reg")
# MasVnrType : NA most likely means no veneer
replace_na(df, "MasVnrType", value="None")
replace_na(df, "MasVnrArea", value=0)
# MiscFeature : data description says NA means "no misc feature"
replace_na(df, "MiscFeature", value="No")
replace_na(df, "MiscVal", value=0)
# OpenPorchSF : NA most likely means no open porch
replace_na(df, "OpenPorchSF", value=0)
# PavedDrive : NA most likely means not paved
replace_na(df, "PavedDrive", value="N")
# PoolQC : data description says NA means "no pool"
replace_na(df, "PoolQC", value="No")
replace_na(df, "PoolArea", value=0)
# SaleCondition : NA most likely means normal sale
replace_na(df, "SaleCondition", value="Normal")
# ScreenPorch : NA most likely means no screen porch
replace_na(df, "ScreenPorch", value=0)
# TotRmsAbvGrd : NA most likely means 0
replace_na(df, "TotRmsAbvGrd", value=0)
# Utilities : NA most likely means all public utilities
replace_na(df, "Utilities", value="AllPub")
# WoodDeckSF : NA most likely means no wood deck
replace_na(df, "WoodDeckSF", value=0)
# CentralAir : NA most likely means No
replace_na(df, "CentralAir", value="N")
# EnclosedPorch : NA most likely means no enclosed porch
replace_na(df, "EnclosedPorch", value=0)
# FireplaceQu : data description says NA means "no fireplace"
replace_na(df, "FireplaceQu", value="No")
replace_na(df, "Fireplaces", value=0)
# SaleCondition : NA most likely means normal sale
replace_na(df, "SaleCondition", value="Normal")
# Electrical : NA most likely means standard circuit & breakers
replace_na(df, "Electrical", value="SBrkr")W dalszej kolejności zamienimy zmienne MSSubClass i
MoSold z numerycznych na kategoryczne, zgodnie z ich
znaczeniem. Zakodujemy także zmienne kategoryczne uporządkowane
(categorical ordinal) z tekstowych na kolejne liczby całkowite.
df = df.replace(
{
"MSSubClass": {
20: "SC20",
30: "SC30",
40: "SC40",
45: "SC45",
50: "SC50",
60: "SC60",
70: "SC70",
75: "SC75",
80: "SC80",
85: "SC85",
90: "SC90",
120: "SC120",
150: "SC150",
160: "SC160",
180: "SC180",
190: "SC190",
},
"MoSold": {
1: "Jan",
2: "Feb",
3: "Mar",
4: "Apr",
5: "May",
6: "Jun",
7: "Jul",
8: "Aug",
9: "Sep",
10: "Oct",
11: "Nov",
12: "Dec",
},
}
)df = df.replace(
{
"Alley": {"None": 0, "Grvl": 1, "Pave": 2},
"BsmtCond": {"No": 0, "Po": 1, "Fa": 2, "TA": 3, "Gd": 4, "Ex": 5},
"BsmtExposure": {"No": 0, "Mn": 1, "Av": 2, "Gd": 3},
"BsmtFinType1": {
"No": 0,
"Unf": 1,
"LwQ": 2,
"Rec": 3,
"BLQ": 4,
"ALQ": 5,
"GLQ": 6,
},
"BsmtFinType2": {
"No": 0,
"Unf": 1,
"LwQ": 2,
"Rec": 3,
"BLQ": 4,
"ALQ": 5,
"GLQ": 6,
},
"BsmtQual": {"No": 0, "Po": 1, "Fa": 2, "TA": 3, "Gd": 4, "Ex": 5},
"ExterCond": {"Po": 1, "Fa": 2, "TA": 3, "Gd": 4, "Ex": 5},
"ExterQual": {"Po": 1, "Fa": 2, "TA": 3, "Gd": 4, "Ex": 5},
"FireplaceQu": {"No": 0, "Po": 1, "Fa": 2, "TA": 3, "Gd": 4, "Ex": 5},
"Functional": {
"Sal": 1,
"Sev": 2,
"Maj2": 3,
"Maj1": 4,
"Mod": 5,
"Min2": 6,
"Min1": 7,
"Typ": 8,
},
"GarageCond": {"No": 0, "Po": 1, "Fa": 2, "TA": 3, "Gd": 4, "Ex": 5},
"GarageQual": {"No": 0, "Po": 1, "Fa": 2, "TA": 3, "Gd": 4, "Ex": 5},
"HeatingQC": {"Po": 1, "Fa": 2, "TA": 3, "Gd": 4, "Ex": 5},
"KitchenQual": {"Po": 1, "Fa": 2, "TA": 3, "Gd": 4, "Ex": 5},
"LandSlope": {"Sev": 1, "Mod": 2, "Gtl": 3},
"LotShape": {"IR3": 1, "IR2": 2, "IR1": 3, "Reg": 4},
"PavedDrive": {"N": 0, "P": 1, "Y": 2},
"PoolQC": {"No": 0, "Fa": 1, "TA": 2, "Gd": 3, "Ex": 4},
"Street": {"Grvl": 1, "Pave": 2},
"Utilities": {"ELO": 1, "NoSeWa": 2, "NoSewr": 3, "AllPub": 4},
}
)W Pandasie zmienne kategoryczne, napisy itp. są typu
object, natomiast numeryczne mają typy z Numpy'a, np.
int64. Oczywiście oba rodzaje zmiennych trzeba przetwarzać
na inne sposoby, więc zapiszemy to od razu. Wyodrębnimy też wektor
docelowy.
Od razu dokonamy też podziału na zbiór treningowy i testowy. Jako że nasz zbiór jest dość mały, to podział będzie 70%-30%.
from sklearn.model_selection import train_test_split
y = df.pop("SalePrice")
categorical_features = df.select_dtypes(include="object").columns
numerical_features = df.select_dtypes(exclude="object").columns
X_train_raw, X_test_raw, y_train, y_test = train_test_split(
df, y, test_size=0.3, random_state=0
)Modele liniowe działają najlepiej, kiedy nasze cechy są odpowiednio przetworzone:
Na start imputujemy wartości brakujące w zmiennych numerycznych i przeskalujemy je do zakresu [0,1]. Potem dokonamy analizy, czy mamy mocno skorelowane zmienne.
Przy przetwarzaniu danych krok po kroku lepiej robić nowe zmienne, żeby dało się odpalić znowu komórkę notebooka. Przy reużywaniu zmiennych mogłaby nastąpić konieczność uruchomienia wszystkiego od początku. Może ci się to przydać w zadaniach.
Pamiętaj, żeby we wszystkich kolejnych zadaniach zastosować taką samą transformację na zbiorze treningowym i testowym. Wszelkie wartości estymujemy na zbiorze treningowym.
Zadanie 2 (0.5 punktu)
Wybierz same zmienne numeryczne. Dokonaj imputacji wartości
brakujących medianą, oraz przeskaluj wartości do zakresu [0,1]. Przydatne klasy:
ColumnTransformer, Pipeline,
SimpleImputer, MinMaxScaler.
Domyślnie ColumnTransformer dopisuje do nazw cech
wyjściowych nazwę transformera jako prefix. Bywa to dość nieczytelne i
problematyczne w dalszym przetwarzaniu danych. Żeby to wyłączyć, przekaż
opcję verbose_feature_names_out=False.
from sklearn.compose import ColumnTransformer
from sklearn.pipeline import Pipeline
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import MinMaxScaler
median_imputer = SimpleImputer(strategy='median')
min_max_scaler = MinMaxScaler()
numerical_transformer = Pipeline(steps=[
('imputer', median_imputer),
('scaler', min_max_scaler)
])
preprocessor = ColumnTransformer(
transformers=[
('num', numerical_transformer, numerical_features)
], verbose_feature_names_out=False
)
X_train_trans = preprocessor.fit_transform(X_train_raw)
X_test_trans = preprocessor.transform(X_test_raw)Sprawdźmy teraz, czy mamy mocno skorelowane zmienne. Typowo przyjmuje się wartość bezwzględną korelacji powyżej progu 0.8/0.9/0.95. Pierwszą weryfikację można zrobić "na oko" na podstawie wykresu typu heatmap. Przekątną oczywiście ignorujemy - każda cecha ma korelację 1.0 z samą sobą.
Jeżeli zauważymy grupy skorelowanych cech, to z każdej grupy trzeba
zostawić jedną cechę. Biblioteka Feature-engine, komplatybilna ze
Scikie-learn, oferuje tutaj wygodną klasę
SmartCorrelatedSelection, która pozwala z każdej grupy cech
zostawić tę o największej wariancji.
Czemu o największej wariancji? Bo taka cecha jest mocno zmienna, a więc niesie jakąś informację. Cechy o bardzo niskiej wariancji praktycznie się nie zmieniają, więc prawdopodobnie nie są zbyt informatywne.
Zadanie 3 (0.5 punktu)
Oblicz korelacje między cechami w zbiorze treningowym, oraz narysuj heatmapę korelacji. Pandas potrafi obliczyć korelacje, a Seaborn - narysować wygodnie heatmapę.
Następnie wykorzystaj Feature-engine do eliminacji mocno skorelowanych cech, z progiem 0.8 (mamy sporo cech, więc eliminujmy, czemu nie), pozostawiając cechy o największej wariancji. Wypisz nazwy wyeliminowanych cech.
correlation_matrix = X_train_trans.corr()
plt.figure(figsize=(20, 16))
sns.heatmap(correlation_matrix, cmap="turbo")
plt.title("Features correlations heatmap")
plt.show()
from feature_engine.selection import SmartCorrelatedSelection
smart_correlated_selector = SmartCorrelatedSelection(
threshold=0.8,
selection_method="variance"
)
smart_correlated_selector.fit(X_train_trans, y_train)
X_train_num = smart_correlated_selector.transform(X_train_trans)
X_test_num = smart_correlated_selector.transform(X_test_trans)
eliminated_features = smart_correlated_selector.features_to_drop_
print("Eliminated features: ", eliminated_features)Eliminated features: ['BsmtFinSF2', 'TotRmsAbvGrd', 'Fireplaces', 'GarageArea', 'GarageCond', 'PoolQC']
Przejdźmy teraz do zmiennych kategorycznych. Jeżeli pewne kategorie występują bardzo rzadko, to prawdopodobnie nie ma sensu ich rozróżniać, i można by je połączyć w jedną, np. "Other" albo "Rare". Eliminacja rzadkich kategorii bardzo zmniejszy wymiarowość po one-hot encodingu, redukując złożoność obliczeniową i pamięciową. Często robi się to nawet głównie z tego powodu, nawet jeżeli lekko pogorszy to wyniki (aczkolwiek często daje nieco lepsze).
Po tym kroku możemy już wykonać one-hot encoding i uzyskać nasze zmienne kategoryczne.
Zadanie 4 (0.5 punktu)
Stwórz ColumnTransformer, wykonujący w pipeline'ie dla
zmiennych kategorycznych (analogicznie jak w zadaniu 1 dla zmiennych
numerycznych):
Dla uniknięcia dummy variable trap, oraz zmniejszenia liczby
wynikowych cech, podaj drop="first" w one-hot
encodingu.
Przyda się RareLabelEncoder z Feature-engine. Żeby
zadziałał dla każdej cechy, podaj mu n_categories=0.
from feature_engine.encoding import RareLabelEncoder
from sklearn.preprocessing import OneHotEncoder
labels_encoder = RareLabelEncoder(n_categories=0,
tol=0.01,
replace_with='Rare')
one_hot_encoder = OneHotEncoder(drop='first', sparse_output=False)
categorical_pipeline = Pipeline(steps=[
('rare_encoder', labels_encoder),
('onehot', one_hot_encoder)
])
categorical_transformer = ColumnTransformer(
transformers=[
('cat', categorical_pipeline, categorical_features)
]
)
X_train_cat = categorical_transformer.fit_transform(X_train_raw)
X_test_cat = categorical_transformer.transform(X_test_raw)Teraz nasze cechy są w pełni przeprocesowane i gotowe do obliczeń. Trzeba je tylko połączyć - zmodyfikuj kod poniżej w zależności od nazw twoich zmiennych.
X_train = pd.concat([X_train_num, X_train_cat], axis=1)
X_test = pd.concat([X_test_num, X_test_cat], axis=1)X_train.shape(2045, 175)
Naszym baseline'owym modelem liniowym, do którego będziemy odnosić bardziej złożone modele, będzie zwykła regresja liniowa (linear regression)
Naszymi metrykami będą zarówno pierwiastek błędu średniokwadratowego (root mean squared error, RMSE), jak i średni błąd bezwzględny (mean absolute error). Oba są w tej samej jednostce, co oryginalna zmienna zależna, ale skupiają się na innych rzeczach:
W praktyce często interesuje nas kilka metryk, a mamy dużo modeli do sprawdzanie. Żeby nie kopiować kodu, warto zrobić funkcję, obliczającą i wypisującą metryki dla danego problemu. Dodatkowo tutaj musimy odwrócić logarytm, który zastosowaliśmy na początku, i takie odwrotne przekształcenia też dobrze jest zawrzeć w takiej funkcji.
from sklearn.metrics import root_mean_squared_error, mean_absolute_error
def evaluate_model(y_train: np.ndarray, y_test: np.ndarray, y_pred_train: np.ndarray, y_pred_test: np.ndarray) -> None:
y_train = np.expm1(y_train)
y_test = np.expm1(y_test)
y_pred_train = np.expm1(y_pred_train)
y_pred_test = np.expm1(y_pred_test)
rmse_train = root_mean_squared_error(y_train, y_pred_train)
mae_train = mean_absolute_error(y_train, y_pred_train)
rmse_test = root_mean_squared_error(y_test, y_pred_test)
mae_test = mean_absolute_error(y_test, y_pred_test)
print(f"Train RMSE {rmse_train:.2f}, MAE: {mae_train:.2f}")
print(f"Test RMSE {rmse_test:.2f}, MAE: {mae_test:.2f}")Jesteśmy teraz gotowi na trening naszych modeli. Na start wytrenujmy zwykłą regresję liniową.
from sklearn.linear_model import LinearRegression
reg_linear = LinearRegression()
reg_linear.fit(X_train, y_train)
y_pred_train = reg_linear.predict(X_train)
y_pred_test = reg_linear.predict(X_test)
evaluate_model(y_train, y_test, y_pred_train, y_pred_test)Train RMSE 17643.27, MAE: 12162.87
Test RMSE 18695.24, MAE: 13107.76
Całkiem nieźle. Dla porównania, bez usuwania skorelowanych cech i grupowania rzadkich cech, RMSE na zbiorze testowym to ponad 21 tysięcy, a model wyraźnie przeucza. Warto wykonywać odpowiedni preprocessing danych!
Mamy tu wyraźną różnicę między RMSE i MAE. Można się temu przyjrzeć i przeanalizować rozkład rezyduów (residuals), czyli błędów y − ŷ.
Zadanie 5 (0.5 punktu)
Narysuj rozkład błędów popełnianych przez model regresji liniowej.
Pamiętaj o odwróceniu logarytmu eksponentą (np.expm1).
Dobierz liczbę binów histogramu tak, żeby wykres był odpowiednio
czytelny. Zaznacz liniami pionowymi średnią oraz medianę. Pamiętaj o
opisaniu osi Y, legendzie i tytule.
Skomentuj rozkład. Jest symetryczny? Może częściej nie doceniamy wartości domu, albo ją przeceniamy, albo może są tu jakiejś pojedyncze outliery? Czy na podstawie tych rezyduów można wytłumaczyć różnicę między MAE i RMSE?
import numpy as np
import matplotlib.pyplot as plt
y_train_orig = np.expm1(y_train)
y_pred_train_orig = np.expm1(y_pred_train)
residuals_train = y_pred_train_orig - y_train_orig
y_test_orig = np.expm1(y_test)
y_pred_test_orig = np.expm1(y_pred_test)
residuals_test = y_pred_test_orig - y_test_orig
plt.figure(figsize=(18, 12))
plt.hist(residuals_train, bins=60, alpha=0.5, label='Train Residuals')
plt.hist(residuals_test, bins=60, alpha=0.5, label='Test Residuals')
plt.axvline(np.mean(residuals_train), color='r', linestyle='--', label='Mean')
plt.axvline(np.median(residuals_train), color='g', linestyle='--', label='Median')
plt.xlabel('Residuals')
plt.ylabel('Frequency')
plt.legend()
plt.title('Distribution of Linear Regression Residuals')
plt.show()
Rozkład błędów popełnianych przez model regresji liniowej wydaje się w miarę symetryczny. Istnieją pojedyncze outliery, które wpływają na zaburzenie symetrii rozkładu.
Różnica między MAE, a RMSE wynika z tego, jak każda z tych metryk traktuje błędy. W przypadku RMSE powodują one większe zaburzenie wyniku.
Mamy pewien overfitting - co prawda niewielki, ale zawsze jest coś do wyeliminowania. Dodatkowo mamy aż 175 cech - jeżeli niektóre z nich mają bardzo niską moc predykcyjną, to redukcja ich wag (lub ich usunięcie) może poprawić nawet i wynik treningowy, i testowy, bo zredukuje szum.
Warto każdy model zapisywać w zmiennej o innej nazwie, żeby później
móc dokonać ich porównania. Przykładowo, reg_linear czy
reg_ridge jest lepszym pomysłem niż samo
reg.
Typowo siłę regularyzacji wybiera się tak, że najpierw szacujemy rząd
wielkości, sprawdzając np. siatkę
[1e-3, 1e-2, 1e-1, 1, 1e1, 1e2, 1e3], a później
zagęszczając ją w okolicy znalezionej optymalnej wartości. Można to
oczywiście robić w kilku krokach.
Zadanie 6 (1 punkt)
Wytrenuj modele z regularyzacją L2 (ridge) oraz L1 (LASSO). Dobierz
siłę regularyzacji, najpierw estymując rząd wielkości, a potem
zagęszczając siatkę. Dla regresji ridge wybieraj model o optymalnym MAE
(argument scoring). Zmierz czas tuningu podczas drugiego
kroku, po zagęszczeniu siatki.
Wypisz optymalne hiperparametry dla obu modeli, oraz ich metryki jakości.
Przydadzą się np.linspace(), RidgeCV oraz
LassoCV.
from sklearn.linear_model import RidgeCV, LassoCV
from time import time
def find_best_alpha(ridge_alphas, lasso_alphas, timer = False):
reg_ridge = RidgeCV(alphas=ridge_alphas, scoring='neg_mean_absolute_error')
reg_lasso = LassoCV(alphas=lasso_alphas)
models = [reg_ridge, reg_lasso]
best_alpha = [None, None]
for i, model in enumerate(models):
start_time = time()
model.fit(X_train, y_train)
if timer:
end_time = time()
print(f"{model.__class__.__name__} tuning time: {round(end_time - start_time, 2)} seconds")
print(f"{model.__class__.__name__} model alpha: {model.alpha_}")
best_alpha[i] = model.alpha_
y_pred_train = model.predict(X_train)
y_pred_test = model.predict(X_test)
evaluate_model(y_train, y_test, y_pred_train, y_pred_test)
print()
print()
return best_alpha, reg_ridge, reg_lasso
reg_values = np.logspace(-3,3,7)
alpha, reg_ridge, reg_lasso = find_best_alpha(reg_values, reg_values)
new_alphas = [np.linspace(alpha[i]/10, alpha[i]*10, 1000) for i in range(2)]
_ , reg_ridge, reg_lasso = find_best_alpha(new_alphas[0], new_alphas[1], timer = True)RidgeCV model alpha: 1.0
Train RMSE 17577.70, MAE: 12172.97
Test RMSE 18481.38, MAE: 12973.79
LassoCV model alpha: 0.001
Train RMSE 20229.58, MAE: 13726.79
Test RMSE 20167.56, MAE: 13954.75
RidgeCV tuning time: 5.98 seconds
RidgeCV model alpha: 1.4081081081081084
Train RMSE 17617.65, MAE: 12196.75
Test RMSE 18468.70, MAE: 12955.05
LassoCV tuning time: 4.1 seconds
LassoCV model alpha: 0.0001990990990990991
Train RMSE 17938.56, MAE: 12366.45
Test RMSE 18500.57, MAE: 12942.30
Dużą zaletą modeli liniowych jest możliwość analizy ważności cech oraz ich kierunku wpływu. W przypadku regresji z regularyzacją jest to szczególnie ciekawe.
Jeżeli cecha "przetrwa" z dużą wagą w ridge regression, to znaczy, że musi być naprawdę ważna, bo inaczej opłacałoby się zmniejszyć jej wagę, żeby zmniejszyć koszt z regularyzacji. Z drugiej strony, w regresji LASSO jeżeli cecha została wyeliminowana, to musiała mieć na tyle małą moc predykcyjną, że modelowi bardziej opłacało się ją usunąć i zmniejszyć koszt z regularyzacji.
Zadanie 7 (0.5 punktu)
W tym zadaniu skorzystaj z optymalnych modeli, wyznaczonych w poprzednim zadaniu.
np.isclose() jej waga jest bliska zeru.Mogą się przydać atrybuty modelu .coef_ oraz
.feature_names_in_. Wszystkie operacje można łatwo wykonać,
konwertując dane na pd.Series.
T = sorted(enumerate(reg_ridge.coef_), key = lambda x: abs(x[1]), reverse=True)
x = [reg_ridge.feature_names_in_[i] for i,_ in T[:20]]
y = [j for _,j in T[:20]]
plt.bar(x, y)
plt.xticks(rotation=90)
plt.title('20 most important features by Ridge Regression')
plt.xlabel('Weight value')
plt.ylabel('Feature name')
plt.show()
lasso_coef = pd.Series(reg_lasso.coef_, index=X_train.columns)
eliminated_features = lasso_coef[np.isclose(lasso_coef, 0)].index.sort_values()
print("Features eliminated by Lasso Regression:")
print(eliminated_features)
Features eliminated by Lasso Regression:
Index(['Alley', 'BsmtCond', 'BsmtHalfBath', 'BsmtUnfSF', 'KitchenAbvGr',
'LowQualFinSF', 'MiscVal', 'PoolArea', 'Street', 'Utilities',
'X2ndFlrSF', 'X3SsnPorch', 'cat__BldgType_2fmCon',
'cat__Condition1_Feedr', 'cat__Condition1_RRAe', 'cat__Condition1_RRAn',
'cat__Condition2_Rare', 'cat__Electrical_Rare',
'cat__Exterior1st_CemntBd', 'cat__Exterior1st_HdBoard',
'cat__Exterior1st_MetalSd', 'cat__Exterior1st_Plywood',
'cat__Exterior1st_VinylSd', 'cat__Exterior1st_WdShing',
'cat__Exterior2nd_VinylSd', 'cat__Exterior2nd_Wd Shng', 'cat__Fence_No',
'cat__Fence_Rare', 'cat__Foundation_Rare', 'cat__Foundation_Slab',
'cat__GarageFinish_No', 'cat__GarageType_Basment',
'cat__GarageType_BuiltIn', 'cat__GarageType_No',
'cat__HouseStyle_1Story', 'cat__HouseStyle_2Story',
'cat__HouseStyle_SFoyer', 'cat__LandContour_Lvl', 'cat__LotConfig_Rare',
'cat__MSSubClass_SC120', 'cat__MSSubClass_SC190',
'cat__MSSubClass_SC60', 'cat__MSSubClass_SC80', 'cat__MiscFeature_Rare',
'cat__MiscFeature_Shed', 'cat__MoSold_Dec', 'cat__MoSold_Mar',
'cat__MoSold_Oct', 'cat__Neighborhood_BrDale',
'cat__Neighborhood_Mitchel', 'cat__Neighborhood_NAmes',
'cat__Neighborhood_SWISU', 'cat__Neighborhood_Sawyer',
'cat__Neighborhood_Timber', 'cat__SaleType_WD '],
dtype='object')
Powyższe algorytmy nie brały pod uwagę skośnego rozkładu zmiennej zależnej, oraz realnego zapotrzebowania na medianę ceny domu. Wypróbujmy zatem algorytmy, które powinny być bardziej odporne na outliery.
Na start wypróbujemy Least Absolute Deviations (LAD) regression. Jako
że jest to przypadek szczególny regresji kwantylowej, to w Scikit-learn
używa się do tego klasy QuantileRegressor z podaniem
odpowiedniego kwantyla.
from sklearn.linear_model import QuantileRegressor
reg_lad = QuantileRegressor(quantile=0.5)
reg_lad.fit(X_train, y_train)
y_pred_train = reg_lad.predict(X_train)
y_pred_test = reg_lad.predict(X_test)
evaluate_model(y_train, y_test, y_pred_train, y_pred_test)Train RMSE 80520.69, MAE: 55705.73
Test RMSE 80920.80, MAE: 55734.18
Błąd jest wręcz tragiczny! W takich przypadkach zawsze na start trzeba się przyjrzeć predykcjom.
y_pred_test[:30]array([11.99535779, 11.99535779, 11.99535779, 11.99535779, 11.99535779,
11.99535779, 11.99535779, 11.99535779, 11.99535779, 11.99535779,
11.99535779, 11.99535779, 11.99535779, 11.99535779, 11.99535779,
11.99535779, 11.99535779, 11.99535779, 11.99535779, 11.99535779,
11.99535779, 11.99535779, 11.99535779, 11.99535779, 11.99535779,
11.99535779, 11.99535779, 11.99535779, 11.99535779, 11.99535779])
Predykcje stałe - nic dziwnego, że błąd jest duży. Scikit-learn domyślnie stosuje lekką regularyzację dla wielu modeli, ale wartość ta może być zbyt duża.
Zadanie 8 (1 punkt)
Dokonaj tuningu siły regularyzacji dla LAD regression, sprawdzając
100 wartości z zakresu od 0 do 1e-3. Wybierz model o
najlepszym MAE z pomocą 5-fold CV. Zmierz czas tuningu, wypisz
znalezione hiperparametry oraz sprawdź jakość finalnego modelu.
Skomentuj różnicę czasową i jakościową względem modeli ridge i LASSO. Czy w tym przypadku warto wykorzystać regresję LAD?
Przyda się GridSearchCV oraz parametr
n_jobs.
from sklearn.model_selection import GridSearchCV
from sklearn.metrics import make_scorer
alpha_values = np.linspace(0, 1e-3, 100)
reg_lad = QuantileRegressor()
param_grid = {
"quantile" : [0.5],
"alpha" : alpha_values}
mae_scorer = make_scorer(mean_absolute_error, greater_is_better=False)
start_time = time()
grid_search = GridSearchCV(reg_lad, param_grid, scoring=mae_scorer, cv=5, n_jobs=-1)
grid_search.fit(X_train, y_train)
end_time = time()
print(f"LAD model alpha: {grid_search.best_params_['alpha']}")
print(f"LAD tuning time: {round(end_time - start_time, 2)} seconds")
y_pred_train_lad = grid_search.predict(X_train)
y_pred_test_lad = grid_search.predict(X_test)
evaluate_model(y_train, y_test, y_pred_train_lad, y_pred_test_lad)LAD model alpha: 0.00047474747474747476
LAD tuning time: 642.73 seconds
Train RMSE 18170.97, MAE: 11649.58
Test RMSE 18189.07, MAE: 12371.05
Wyniki LAD Regression są zbliżone do Ridge i LASSO, natomiast czas wykonania dużo dłuższy z powodu konieczności obliczenia układów równań programowaniem liniowym. Zatem w tym przypadku nie warto korzystać z LAD Regression.
Zamiast używać regresji LAD, opartej o programowanie liniowe, można też wykorzystać Huber regression. Jest ono mało czuła zarówno na ϵ (punkt przejścia z MSE do MAE), jak i α (siłę regularyzacji L2), więc często nie wymaga tuningu hiperparametrów.
Zadanie 9 (0.5 punktu)
Wytrenuj model Huber regression. Sprawdź jakość modelu. Skomentuj różnicę czasową i jakościową względem LAD regression (tutaj oczywiście bez tuningu). Z którego modelu skorzystałbyś w praktyce?
Domyślna liczba iteracji solwera, 100, jest zbyt mała dla naszej liczby cech, i będą z tego wynikać ostrzeżenia o braku zbieżności, więc zwiększ ją do 2000.
from sklearn.linear_model import HuberRegressor
reg_huber = HuberRegressor(max_iter=2000)
start_time = time()
reg_huber.fit(X_train, y_train)
end_time = time()
print(f"Huber tuning time: {round(end_time - start_time, 2)} seconds")
y_pred_train_huber = reg_huber.predict(X_train)
y_pred_test_huber = reg_huber.predict(X_test)
evaluate_model(y_train, y_test, y_pred_train_huber, y_pred_test_huber)Huber tuning time: 2.34 seconds
Train RMSE 18582.51, MAE: 11701.59
Test RMSE 17994.72, MAE: 12252.40
Huber Regression daje porównywalne wyniki (może nawet trochę lepsze), co Ridge i LASSO w stosunkowo podobnym czasie. Na plus dodatkowo jest to, że nie zwraca tak uwagi na outliery.

Na rynku działa firma oferująca serwis do sprzedaży nieruchomości. CEO był na prezentacji o AI, zaczął używać ChatGPT i uznał, że też musicie mieć taki potężny ML. Zostałeś zatrudniony jako (pierwszy i jedyny) data scientist, zebrałeś historyczne dane, i nawet wdrożyłeś pierwszy model regresji liniowej do estymacji wartości nieruchomości.
Pewnego dnia Product Owner przychodzi do ciebie z listą "absolutnie krytycznych" wymagań (propozycja nie do odrzucenia):
Zadanie 9 (4 punkty)
ColumnTransformer. Ma do niego wejść cale
X_train_raw, a wyjść gotowy output.ames_transformer.joblib,
price_estimator.joblib,
low_price_estimator.joblib,
high_price_estimator.joblib.Poniżej przygotowano funkcję testową, która ładuje transformer i estymatory, oraz oblicza predykcje dla trzech przykładowych domów ze zbioru testowego (taniego, przeciętnego i drogiego).
Skomentuj - czy twoim zdaniem finalne modele, podające przeciętną wartość i widełki, są subiektywnie dobrej jakości?
from feature_engine.outliers import Winsorizer
import joblib
capper = Winsorizer(capping_method='gaussian',
tail = "right",
fold=3)
num_overall_pipeline = Pipeline([("imput", median_imputer),
("scaling", min_max_scaler),
("winsorization", capper),
("selection", smart_correlated_selector)])
category_pipeline = Pipeline([("rare_labels", labels_encoder),
("one_hot_encoder", one_hot_encoder)])
overall_transformer = ColumnTransformer([("numerical transform", num_overall_pipeline, numerical_features),
("category transform", category_pipeline, categorical_features)],
verbose_feature_names_out=False)
X_train_overall = overall_transformer.fit_transform(X_train_raw)
X_test_overall = overall_transformer.transform(X_test_raw)
joblib.dump(overall_transformer, "ames_transformer.joblib")['ames_transformer.joblib']
huber_true_price = HuberRegressor(max_iter=2000)
huber_true_price.fit(X_train_overall, y_train)
huber_low_price = HuberRegressor(max_iter=2000)
huber_low_price.fit(X_train_overall, np.log1p( np.expm1(y_train) * 0.35))
huber_high_price = HuberRegressor(max_iter=2000)
huber_high_price.fit(X_train_overall, np.log1p( np.expm1(y_train) * 0.65))
joblib.dump(huber_true_price, "price_estimator.joblib")
joblib.dump(huber_low_price, "low_price_estimator.joblib")
joblib.dump(huber_high_price, "high_price_estimator.joblib")['high_price_estimator.joblib']
import joblib
def test_models() -> None:
transformer = joblib.load("ames_transformer.joblib")
price_reg = joblib.load("price_estimator.joblib")
low_price_reg = joblib.load("low_price_estimator.joblib")
high_price_reg = joblib.load("high_price_estimator.joblib")
# select houses from test set around 25th, 50th and 75th percentile
prices = np.expm1(y_test).astype(int)
for price_type, quantile in [("Cheap", 0.25), ("Average", 0.5), ("Expensive", 0.75)]:
idx = prices[prices == prices.quantile(quantile, interpolation="nearest")].index[0]
X = X_test_raw.loc[idx].to_frame().T
y = prices[idx]
float_to_int = lambda x: int(np.round(np.expm1(x.item())))
X_transformed = transformer.transform(X)
y_pred = float_to_int(price_reg.predict(X_transformed))
y_low_pred = float_to_int(low_price_reg.predict(X_transformed))
y_high_pred = float_to_int(high_price_reg.predict(X_transformed))
print(f"{price_type} house")
print(f"True price: {y}$")
print(f"Estimated price: {y_pred}")
print(f"Price brackets: {y_low_pred} - {y_high_pred}")
print("Features:")
with pd.option_context("display.max_columns", None):
display(X)
test_models()Cheap house
True price: 127999$
Estimated price: 133395
Price brackets: 46691 - 86708
Features:
| MSSubClass | MSZoning | LotFrontage | LotArea | Street | Alley | LotShape | LandContour | Utilities | LotConfig | LandSlope | Neighborhood | Condition1 | Condition2 | BldgType | HouseStyle | OverallQual | OverallCond | YearBuilt | YearRemodAdd | RoofStyle | RoofMatl | Exterior1st | Exterior2nd | MasVnrType | MasVnrArea | ExterQual | ExterCond | Foundation | BsmtQual | BsmtCond | BsmtExposure | BsmtFinType1 | BsmtFinSF1 | BsmtFinType2 | BsmtFinSF2 | BsmtUnfSF | TotalBsmtSF | Heating | HeatingQC | CentralAir | Electrical | X1stFlrSF | X2ndFlrSF | LowQualFinSF | GrLivArea | BsmtFullBath | BsmtHalfBath | FullBath | HalfBath | BedroomAbvGr | KitchenAbvGr | KitchenQual | TotRmsAbvGrd | Functional | Fireplaces | FireplaceQu | GarageType | GarageYrBlt | GarageFinish | GarageCars | GarageArea | GarageQual | GarageCond | PavedDrive | WoodDeckSF | OpenPorchSF | EnclosedPorch | X3SsnPorch | ScreenPorch | PoolArea | PoolQC | Fence | MiscFeature | MiscVal | MoSold | YrSold | SaleType | SaleCondition | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 195 | SC50 | RM | 50.0 | 6000 | 2 | 0 | 4 | Lvl | 4 | Inside | 3 | BrkSide | Norm | Norm | 1Fam | 1.5Fin | 6 | 6 | 1940 | 1950 | Gable | CompShg | MetalSd | MetalSd | None | 0.0 | 3 | 4 | CBlock | 3 | 3 | 0 | 2 | 264.0 | 1 | 0.0 | 308.0 | 572.0 | GasA | 5 | Y | FuseA | 848 | 348 | 0 | 1196 | 0.0 | 1.0 | 1 | 1 | 3 | 1 | 3 | 6 | 8 | 2 | 4 | Detchd | 1973.0 | Unf | 2.0 | 576.0 | 3 | 3 | 2 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | No | No | 0 | Mar | 2010 | WD | Normal |
Average house
True price: 157999$
Estimated price: 151862
Price brackets: 53153 - 98714
Features:
| MSSubClass | MSZoning | LotFrontage | LotArea | Street | Alley | LotShape | LandContour | Utilities | LotConfig | LandSlope | Neighborhood | Condition1 | Condition2 | BldgType | HouseStyle | OverallQual | OverallCond | YearBuilt | YearRemodAdd | RoofStyle | RoofMatl | Exterior1st | Exterior2nd | MasVnrType | MasVnrArea | ExterQual | ExterCond | Foundation | BsmtQual | BsmtCond | BsmtExposure | BsmtFinType1 | BsmtFinSF1 | BsmtFinType2 | BsmtFinSF2 | BsmtUnfSF | TotalBsmtSF | Heating | HeatingQC | CentralAir | Electrical | X1stFlrSF | X2ndFlrSF | LowQualFinSF | GrLivArea | BsmtFullBath | BsmtHalfBath | FullBath | HalfBath | BedroomAbvGr | KitchenAbvGr | KitchenQual | TotRmsAbvGrd | Functional | Fireplaces | FireplaceQu | GarageType | GarageYrBlt | GarageFinish | GarageCars | GarageArea | GarageQual | GarageCond | PavedDrive | WoodDeckSF | OpenPorchSF | EnclosedPorch | X3SsnPorch | ScreenPorch | PoolArea | PoolQC | Fence | MiscFeature | MiscVal | MoSold | YrSold | SaleType | SaleCondition | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 2538 | SC80 | RL | 80.0 | 9600 | 2 | 0 | 4 | Lvl | 4 | Inside | 3 | NAmes | Norm | Norm | 1Fam | SLvl | 5 | 7 | 1967 | 1967 | Gable | CompShg | MetalSd | MetalSd | BrkFace | 140.0 | 3 | 3 | PConc | 3 | 3 | 2 | 5 | 602.0 | 3 | 402.0 | 137.0 | 1141.0 | GasA | 4 | Y | SBrkr | 1141 | 0 | 0 | 1141 | 1.0 | 0.0 | 1 | 0 | 3 | 1 | 3 | 6 | 8 | 0 | 0 | Attchd | 1967.0 | Unf | 1.0 | 568.0 | 3 | 3 | 2 | 0 | 78 | 0 | 0 | 0 | 0 | 0 | No | No | 0 | Jul | 2006 | WD | Normal |
Expensive house
True price: 209499$
Estimated price: 190714
Price brackets: 66749 - 123940
Features:
| MSSubClass | MSZoning | LotFrontage | LotArea | Street | Alley | LotShape | LandContour | Utilities | LotConfig | LandSlope | Neighborhood | Condition1 | Condition2 | BldgType | HouseStyle | OverallQual | OverallCond | YearBuilt | YearRemodAdd | RoofStyle | RoofMatl | Exterior1st | Exterior2nd | MasVnrType | MasVnrArea | ExterQual | ExterCond | Foundation | BsmtQual | BsmtCond | BsmtExposure | BsmtFinType1 | BsmtFinSF1 | BsmtFinType2 | BsmtFinSF2 | BsmtUnfSF | TotalBsmtSF | Heating | HeatingQC | CentralAir | Electrical | X1stFlrSF | X2ndFlrSF | LowQualFinSF | GrLivArea | BsmtFullBath | BsmtHalfBath | FullBath | HalfBath | BedroomAbvGr | KitchenAbvGr | KitchenQual | TotRmsAbvGrd | Functional | Fireplaces | FireplaceQu | GarageType | GarageYrBlt | GarageFinish | GarageCars | GarageArea | GarageQual | GarageCond | PavedDrive | WoodDeckSF | OpenPorchSF | EnclosedPorch | X3SsnPorch | ScreenPorch | PoolArea | PoolQC | Fence | MiscFeature | MiscVal | MoSold | YrSold | SaleType | SaleCondition | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 192 | SC75 | RL | 0.0 | 7793 | 2 | 0 | 3 | Bnk | 4 | Corner | 3 | BrkSide | Norm | Norm | 1Fam | 2.5Unf | 7 | 7 | 1922 | 2005 | Gable | CompShg | Wd Sdng | Wd Sdng | None | 0.0 | 3 | 3 | BrkTil | 4 | 3 | 0 | 4 | 474.0 | 1 | 0.0 | 634.0 | 1108.0 | GasA | 3 | N | FuseA | 1160 | 908 | 0 | 2068 | 0.0 | 0.0 | 1 | 1 | 3 | 1 | 4 | 8 | 8 | 1 | 4 | Detchd | 1928.0 | Unf | 1.0 | 315.0 | 3 | 3 | 2 | 0 | 0 | 60 | 0 | 0 | 0 | 0 | No | No | 0 | May | 2010 | WD | Normal |
Wybrałem Huber Regression, ponieważ nie zwraca zbyt dużej uwagi na outliery, a modele czasami niedoceniają drogich domów lub przeceniają tanie.
Product Owner przychodzi do ciebie ponownie. Klienci skarżą się, że
nie rozumieją, czemu model przewiduje taką, a nie inną wartość dla ich
domu. Przykładowo, jeżeli to ogólny stan domu (OverallCond)
czy garażu (GarageCond) obniżają cenę, to byliby gotowi je
wyremontować dla lepszego zysku.
Jest to problem wyjaśnialnego AI (Explainable AI, XAI). O ile cały model w przypadku regresji liniowej jest bardzo transparentny dzięki wagom cech i liniowym predykcjom, to już zrozumienie, czemu dana predykcja była taka, a nie inna, to już problem lokalnej wyjaśnialności (local explainability), gdzie tłumaczymy pojedynczą predykcję, dla danego zestawu cech.
Najpopularniejszym algorytmem jest tutaj SHapley Additive Explanations (SHAP). Opiera się on o tzw. wartości Shapleya, mające bardzo solidne podstawy w kooperatywnej teorii gier. Opiera się na idei, że wartości cech grają w grę: jedne obniżają predykcję, a drugie podwyższają, a finalna predykcja modelu to ich konsekwencja. Każda cecha to gracz, ale jedna wnosi więcej, a inna mniej. Idea algorytmu SHAP to uznanie, że "wartościowi gracze" to wpływowe cechy, które mocno przyczyniły się do danej decyzji modelu.
Obliczanie wartości Shapleya w ogólnym przypadku ma złożoność eksponencjalną. Algorytm SHAP zaproponował efektywne ich przybliżanie dla dowolnych algorytmów ML. Dla niektórych klas modeli, przede wszystkim liniowych oraz drzewiastych, istnieją szczególnie efektywne metody, które w czasie wielomianowym obliczają dokładne wartości Shapleya. Wyjaśnienia, jak działa SHAP: link 1, link 2, link 3.
Z pomocą biblioteki shap oblicz wartości Shapleya dla
przykładowych domów z ostatniego zadania i wyświetl je na tzw. force
plot. Przyda ci się tutaj LinearExplainer (dokumentacja)
oraz plots.force() (dokumentacja).
Porównaj ze sobą dwie metody estymacji: "interventional" oraz
"correlation_dependent" (opisane w dokumentacji). Możesz je
zaimplementować za pomocą maskers.Independent oraz
maskers.Impute. Jako zbiór do porównania (tzw. background
samples) użyj całego zbioru treningowego.
Skomentuj, czy wyjaśnienia mają twoim zdaniem sens i czy są intuicyjne. Która metoda estymacji daje subiektywnie lepsze wyniki?
import shap
from shap import maskers
transformer = joblib.load("ames_transformer.joblib")
price_reg = joblib.load("price_estimator.joblib")
background = X_train_overall
masker_independent = maskers.Independent(background)
masker_impute = maskers.Impute(background)
explainer_independent = shap.LinearExplainer(price_reg, masker_independent)
explainer_impute = shap.LinearExplainer(price_reg, masker_impute){"model_id":"10cd15c0ac244efe9428676a07414c3c","version_major":2,"version_minor":0}shap.initjs()
shap.plots.force(explainer_independent.expected_value, shap_values_independent[0])shap.initjs()
shap.plots.force(explainer_impute.expected_value, shap_values_impute[0])shap.initjs()
shap.plots.force(explainer_independent.expected_value, shap_values_independent)shap.initjs()
shap.plots.force(explainer_impute.expected_value, shap_values_impute)Wyjaśnienia wydają się być dość intuicyjne. Metoda Correlation Dependent wydaje się być lepsza w tym przypadku, ponieważ zakłada korelację cech i uwzględnia w obliczeniach. Metoda Interventional zakłada ortogonalność cech przez co może przekłamywać wyniki.